flASSACHUSETTS INSTITUTE OF TECHNOLOGY A.I. LABORATORY Artificial Intelligence October 1973 demo. No. 235 ON LIGHTNESS Berthold K. P. Horn The intensity at a point in ^n image is the product of the reflectance at the corresponding object point and the intensity of illumination at that point. Lie are able to perceive llghtneee, a quantity closely correlated ulth reflectance. Hon then do Me eliminate the component due to illumination from the Image on our retina? The tuo components of image Intensity differ in their epatial distribution* A method Is presented here which takes advantage of thle to compute lightness from image intensity in a layered, parallel fashion. The method is developed for a restricted class of images first used by Land in presenting his retinex theory of color. In thle theory the problem of color perception is reduced to one of judging black and white lightness on three images taken in different parts of ths visual spectrum. The method described h*re fills the need for a lightness judging procees. The theory has implications for potential special purpose hardware In Image sensing devices. It should also bs of interest to cognitive peychologiete elnce it can explain certain effects observed in the human vleu*l system ae well ae predict net* ones. Further, the theory provides neuro-phyelologlsts with suggsstlons about the function of certain etructuree In the primate retina. This report describes research done at the Artificial Intelligence Laboratory of the Massachusetts Institute of Technology. Support for the laboratory's artificial intelligence research is provided in part by the Advanced Research Projects Agency of the Department of Defense under Office of Naval Research contract N000U-70-A- OM2-OM5. Reproduction of thle document. In whole or in part, la permitted for any purpose of the United States Government* 1. Revieu. 3 1.1 Theories of Color Perception. 3 1.1.1 Tri-Stlmulue Theory. * 1.1.2 Color Conversion. 4 1.2 Land's Retinex Theory. S 1.2.1 Lightness Judging. 6 1.2.2 fllni-HorJd of Mondriane. 6 1.3 Uhy Study the One-dimensional Case? 7 1.3.1 Notation. 8 1.4 One-Dimensional Method - Continoue Case. 9 1.4.1 One-Dlmenslonal Continoue Method: Details. 9 1.4.2 Normalization. 12 1.5 One-Dimensional Method - Diecrete Case. 13 1.5.1 One-Dimensional Di ecrete flethod: Details. 13 1.5.2 Selecting the Threshold. IS 1.5.3 Accuracy of the Reconstruction. 17 1.5.4 Generalizations. 17 1.5.5 Physical Models of the One-Dimensional Discrete Process* 18 2. Lightness in Two Dimensional Images. 19 2.1 Tuo-Dineneional Method - Continoue Case. 19 2.1.1 Applying the Laplaclan to a Mondrian. 19 2.1.2 Inverse of the Laplacian Operator, 29 2*1.3 Uhy one can use the Convolutions! Inverse. 23 2.1.4 Normalization. 23 2.1.5 Tuo-D i mens I ona I Continous Methodr Detal Is- ♦ - * 24 2.2 Tuo-Dimeneional Method - Discrete Caee. 25 2.2.1 Teeselatlon of the Image Plane. 25 2.2.2 Diecrete Analogue of the Laplacian. 26 2.2.3 Inverse of the Diecrete Operator. 27 2.2.4 Computational Effort and Simplification. 29 2,2*5 Tuo-Dimeneional Discrete Methods Oetatle. 38 2.2.6 Simplicity of the Inverse. 32 2.2.7 Iterative Methods of Solution, « 32 2.2.8 Convergence of Iterative and Feed-back Schemes. 33 2.2.9 Setting the Threehold. 34 2.2.18 Some Notes on This Msthod. 35 2.2*11 Dynamic Range Reduction. 36 2.2.12 A Frequency Domain Interpretation. 36 2.3 Physical Mode J e. 37 2.3.1 A Diecrete Physical Modsl. 38 2.3.2 A Feed-back Schsme for the Inverse. 48 2.4 Limitations of the Simple Scheme Presented. 42 2.5 Computer Simulation of the Discrete Method. 44 2.5.1 Form of Inveree ueed in the Computer Simulation. 4S 3. Impl icatione and Concluelone* 47 3.1 Parallel Image Proceeeing Hardware. 47 3.2 Cognitive Peychology, *7 3.3 Neuro-phyeiology. 4 * 3.4 Conclusion. ** ON LIGHTNESS 2 Acknoulegaente I wish to thank Dr. David flarr for a numbsr of most Interesting dlacueeione and for encouraging na to commit this msthod to papar. Ha first observed that the primate retina appears to have Juet the right detailed etructure to implement the functions here described. LIGHTNESS* Definition The relative degree to which an object reflecte light. The Random House Dictionary The attribute of object colors by which the object appears to reflect or transmit more or less of the incident light* Uabster's Seventh Nsu Collaglata Qictionary Preview Part 1 Is a review of the relevant Information relating to color vleion and lightness. Thie Includes a discussion of the Land ratine* model in a form suitable for the developments of the next part. In part 2, Land's one-dimensional operation will be extended to two-dimensional Images. The method depends on a layered, parallel computation suggestive of both biological and artificial implementations. In part 3 some of the implications are explored and the information on the new image processing technique is summarized. ON LIGHTNESS 3 1. Review, 1.1 Theories of Color Perception. Thar* has always been great interest in how we perceive colors and numerous explanations have been forwarded (Newton 1704, Goethe) 1818, Young 1820, Maxwell 1856, Helisholt* 1867, Hering 1875). The huaan perceptual apparatus Is remarkably eucceeful in coping with large variations in the i I lueinatlon. The colors us perceive are closely correlated ulth the surface colours of the objecte viewed, deepite large temporal and epatial differences in color and intensity of the Incident light. Thie is eurprieing eince we cannot eenee reflectance directly. The light inteneitg at a point in the iaage le the product of the reflectance at the corresponding object point and the intensity of illumination at that point - aside from a constant factor that depende on the optical arrangment. There must then be soee difference between theee two components of Image intensity which allows ue to dlecount the effect of one. The two componente differ In their spatial distribution. Incident light inteneitg will usually vary smoothly, with no di econt inui t ies, uhlls reflectance will have eharp diecontinul ties at edges where objecte adjoin. The rsflsctance being relatively constant between euch edgea. ON LIGHTNESS k 1*1.1 Trl-Stlmului Theory* Some facto about how we see color are fairly uell established. It appears that us have three kinds of sensors operating In bright 1 1 lull n*tl on, Hith peak aanai tivi tlee in different parts of the visible spectrum. This is why it takes exactly three colors in additive mixture to natch an unknoun color. I4h1 le It la a bit tricky to Measure the sensitivity curvss of the three sensors directly, a linear transform of these curves has been known accurately for some tine (Brlndley 19GB). Theee curves, ceiled the standard observer curves, are sufficient to allow one to predict color matches made by subjecte with normal colour vieion (Hardy 193G1 . The elmplest theory of color perception then amounte to locally comparing the outputs of three such sensors and assigning colour on thle basis (Young 1820, Helmholtz 18671. This however totally falls to explain the observed color contancy. Perceived color does not depend directly on the relative amounts of light msasured by the three sensors (Land 1959, Lettvin 1967). 1.1.2 Color Conversion. A number of attempts have been mads to patch up this theory under the rubrics of "discounting of the lllumlnanf, "contraet effect adjuetment" and "adaptation". The mors complicated theories are based on models with large numbers of paramstsrs which era adjusted according ON LIGHTNESS S to empirical data IHelson 1938 & 1943, Judd 1348 & 1952, Richard* 1971). Theee theories are at leaet partially effective in predicting human color perception when applied to simple arrangements of stimuli similar to those used in determining the parameters. The parameters depend strongly on the data and alight experimental variations will produce large fluctuations in them. This Is a phenomena familiar to numerical analysts fitting curves to data when the number of parameters is large. These theories are lacking in parsimony and convincing physiological counter-parts. Lettvln has demonstrated the hopeleesneae of trying to find flxsd transformations from locally compared output of sensors to perceived color {Lettvln 1967). 1.2 Land*s Retinex Theory. Another theory of color perception is embodied In Land'e retinex model (Land 1953. 1964 & 19711. Land proposes that the three eete of seneore are not connected locally, but instead are treated ae If they represent points on three separate images. Processing is performed on each euch image separately to remove the component of Intensity due to illumination gradient. Such processing is not msrely an added frill but ie indi sponsible to color perception in the face of the variability of the i I lumination. ON LIGHTNESS 6 1*2.1 Lightness Judging* In essence a judge of lightness processes each image* Lightness is the perceptual quantity closely correlated with surface reflectance. Only after this process can the three images be compared to reliably determir>e colors locally. It remains to mechanize this process. It Mould appeal to Intuition If this process could be carried out in a parallel fashion that doss not depend on previous knowledge of the scene viewed* This is because colors are so immediate, and seldom depend on one's interpretation of the scene. Colore will be seen even when the picture makes no ssnse in terns of previous experience. Also, color Is eeen at every point In an Image. 1*2.2 lllni-world of tlondrians. In developing and explaining his theory Land needed to postpone dealing ulth the full complexity of arbitrary scenes* He eelected a particular claes of objscts as Inputs, modelled after the paintings of the turn-of-the-century Dutch artist Pieter Cornells (londrian. These scenes are flat areas divided into sub-regions of uniform matte color* Problems such as those occasionsd by shadows and specular reflection are avoided In this way. One also avoids shading; that is, the variation In reflectance with the orientation of the surface In respect to the sensor and the light-source. For Hondrians* lightness is considered to be a function of reflectance* ON LIGHTNESS 7 llondrlans ara usually made of polygonal regione mth straight aides — for the development here however the edges may be curved. In the world of Mondrians one finds that the reflectance has sharp discontinuities wherever regions meet, being constant inaida each region. The illumination, on tha other hand, varies anoothlg over the image, 1*3 Uhy Study the One-dimensional Case? Imagee MrB tuo-dlmeneional and usually sampled at discrete points. For hlatoric reasons and intuitive simplicity tha results will flret be developed in one dimension, that is with functions of one variable* Similarly, contlnoue functlone uill be used at flret eince they aJIOM a cleaner separation of the two componsnts of image intensity and iiluetrate nore clearly the concepts involved. Use Hill be Bade of analogies between the one-dimensional and two- dimensional casss as well ae the continous and discrete ones. The final procaee discussed for processing image intensities Is tuo-dlmenelonal and dlecrete. A number of physical lmplenentatlons for this scheme ar9 suggested. The proceae hIII be looked at from a number of points of view: partial dl f ferent ial equat ione, I inear systems, four ter transforms and convolutions, difference equations, iterative solutions, feed-back schemes and physical models. ON LIGHTNESS 8 1*3.1 Notation. The following notation mi I I be used. a* Intensity of incident illumination at a point on the object* r' Reflectance at a point of the object* p* Intensity at an image point* Product of •' and r'. a t r. pi Logarithms of e\ r* and p' respectively. d Result of applying forward or differencing operator to p. t Result of applying threshold operator to d. I Result of applying inverse or summing operator to t* Simple derivative operator in one dimension* T Contlnous threshold operator, discards finite part. 1 Simple integration operator in one dimension, L Laplacian operator - eum of second partial derivatives. G' Inveree of the Laplacian, convolution uith (1/2 t ) log c d/r). D*, T*, I*, L* and Gfri Dlacrete analogues of 0, T, I, L and G. The output It will not be called lightness since there le probably not yet a generally acceptable definition of this term. It is however intended to be monotonical ly related to lightness. Note that I le ON LIGHTNESS 3 related to the logarithm of reflectance, uhile the perceptual quantity le perhape more closely related to the square-root of reflectance. 1.4 Qne-Oimensional Plethod - Continous Case. Land invented a simple method for separating the Image components in one dimension* First one takes logarithms to convert the product Into a sum* This Is followed by differentiation. The derivative uill be the sum of the derivatives of the two components. The edges uill produce sharp pulses of arsa proportional to the Intensltg eteps between regione - uhi le the spatial variation of II I urn (nation uill produce only f i nf to values overyuhere. Noh if ono diacarda all finite yatuae, one It left ulth the pulses and hsnce the derivative of lightness* Finally one undoee the differentiation by simple integration. 1.4.1 One-Dimensional Continous Method! Details. Ue have the following. Let r'UI be the reflectance of the object at the point corresponding to the image point x. Let eMx) be the intensity at this objact point. Let p* (x) be their product, that ie, the intensity recorded In the imags at point x. Note that e' (x) and r* (k) %r9 poei tive, p* <x) - e' (x) * rMx) .« t 1 T dlx T tlxj ■ H FIGURE 1: Processing steps in the one-dimensional continous case* ON LIGHTNESS 11 Now let p(x) be the logarithm of p' (x) and bo oni ptx) - six) + r(x) Note that six) is continous and that r(x) has some finite diacont inui ties. Let repcssent differentiation with reapect to x* dfx) - D(p(x>) - 0(a(x)) + D(p(x)) Now, D(s(x)) Hill be finite everywhere, while 0(r(x)) Mill be zero aalde from a number of pulses - which carry ail the information. Each pulae Mill correspond to an edge between regions and have area proportional to the intensity atep. If now one ■thresholds'* and discards all finite parte, one getei Km) - T(D(p(x)>) ■ D(r(»)) To obtain r(x) one only has to invert the differentiation, that is, integrate. Let I represent integration with respect to x f then (I) - D andi IlKl - I(T(D(p(x))) - rtx) + c One can give a convolutions! Interpretation to the above, since differentiation corresponds to convolution with a pulse-pair, one ON LIGHTNESS 12 negative and one positive, each of unit area. Integration corresponds to convolution with the unit step function. 1.4.2 Normalization* The reeult is not unique because of the constant introduced by the integration. The ztro (spatial) frequency term has been loet in the differentiation, so cannot be reconstructed. Thle le related to the fact that one does not knoM the overall level of Illumination and hence cannot tell whether en object appears dark because it is grey or because the level of illumination is loti. One can normalize the result If one assumes that there are no light •oyrces in the field of view and no floureecent colore or specular reflections, Thle is certainly the case for the flondrians, Perhape the beet way of normalising the result Is to simply assume that the highest value of lightness corresponds to white* or total reflectance in the Lambertian sense. This normalization will lead one astray if the image doee not contain a region corresponding to a white patch in the ecene, but this is the beet one can do. Other normalization techniques might involve adjusting weighted local averages, but thle would then no longer amount to reconstruction of reflectance. ON LIGHTNESS 13 l.S Qne-Oimenaional Hethod - Discrete Case. So far ue have assumed that tha image intensity uan a continoue function. In retinas found In animal a or artificial onss constructed out of dlacrata componsnte. iniagaa ara only saitplsd at discrete points'. So ona haa to find diacrats analogues for the operations we have been using. Perhaps tha simplest ar* first differences and summation as analogues of differentiation and integration respectively. This is not to say that other approxi mat ions could not be used equally wall. To use the new operators one goes through essentially tha sane process as before, except that now all values in the differenced imaga ara finite. This has the effect of forcing one to choose a threshold for the thresholding function. Both components of image intensity produce finite values after tha differencing operation. The component due to the edges in the reflectance is hopefully quits large compared to that due to illumination gradient. One has to find a level that Mill suppress the illumination gradient inside regions, uhi la permitting the effecte due to edges to reaain, 1*8.1 One-Dimensional Discrete flethod: Details. Let rl be the reflectance of the object at the point corresponding to the image point 1. Let s\ be the incident light intensity at this object point. Let p\ be their product* that la tha lntenelty In tha image at point I. - .* D j=l ii «; FIGURE 2: Processing steps In the one-dimensional discrete case, ON LIGHTNESS IS Now 1st pj b* th* log of p- and bo on. Lot Oft and Ift 0* !"• operators corresponding to taking first differences and summation respect Ively. Note that (Ift)"' - D*. Pi- •l+n d i- p*h- p; (d - d * ( p ,j t ; - d; If !d,[ < «. tlio -£' ll- > t K - K-0 1*5.2 Selecting the Threehotd. Uhat determines the threshold? It mutt be smaller than the emalleet inteneitg etep between regions. It ituet on the other hand be larger than values produced by first differencing the maximum Illumination gradients. Real images are noisy and the threshold ehould be large enough to eliminate this noise inside regions* The spacing of the seneor cells must also be taken Into account* GN LIGHTNESS 16 As this spacing becossv smaller, the contribution due to Illumination gradients decreases, uhlls the coaponsnt due to the edges remains constant. A limit is reached uhen the component due to illumination gradients falls below that due to noise or when the optical properties of the Imaging system begin to have a deliterlous effect. In all imaging systems an edge Is spread over a finite distance due to diffraction and uncorrected abberatione. The spacing of seneore should not be much smaller than this distance to avoid reducing the component due to edges In the differenced iaags. Let u be the radius of the point-spread-function of the optical system and h the spacing of ths ssnsor celle. Let g* be the smallest atop in the logarithm of raflectanca In tha acana. Than define tha affective minimus step asi g - g* a mind, h/2u) Let a be the largest slope due to illumination gradient and * the root- mean-equare noise-amplitude. The noise will exceed a value 3 only «3% of the time. Choose the threshold e as follows* a < g • > a * h e > 3\/2 a ON LIGHTNESS 17 1.6.3 Accuracy of the Reconstruction. In thA contlnoue case one can exactly reconstruct tha reflectance, aeide from a constant. Us are not eo fortunate here, even If ue Mlect a threshold according to the above criteria. This la because the valuee at the edges contain email contributions due to Illumination gradient and nolee. A alight Inaccuracy In the reconetructlon ulll reeult. Thle error le minimized by making the esneor cell spacing very fine, optimally of a elze comaeneurate with the optical reeolutlon of the device. The effect of nolee can also be minimized by integrating over tlM. Note that the reconstruction le aore accurate when there are fen edgee, elnce It le at the edges that the error effects appear. Ulth many edgee the Illumination gradient begins to "show through - . 1*5.4 General izat lone. So far we have dealt ulth constant sensor epaclng. Clearly aa long ae the aame epaclng le ueed for both the differencing and the summing, the cell epaclng can be arbitrary and has little effect on the reconetruction elnce It does not enter Into the equatlone. ■ Similarly ue have chosen first differences as the dlecrete analogue for differentiation. Ue could have chosen some other uelghted difference and developed • eultable Inverse for It. Thle Inveree of course uould no longer be eummatlon but cin be readily obtained using ON LIGHTNESS 18 techniques developed for dealing with difference equetlone (Rlchtmeyer 1957, Garabedlen 19G4i. 1*5.5 Physical rtodels of the Ons-Oimensional Discrete Proceee. One can invent a number of physical modele of the above operatlone. A simple reeistive netuork ulll do for the summation proceee for example. Land has implemented a small circular "retina" ulth about 1G eensora. This model employs electronic componente to perform the operations of taking logarithms, differencing, thresholding and summing. Land hee tried to extend his one-dimensional method to Images* by covering the image ulth paths produced by a random ualk procedure and applying methods like the above to each of these paths. Uhllo this produces reeulte. It esens unsatisfactory from the point of view of euggeetlng possible neuro-physlologlcal structureet neither does It lend itself to efficient Implementation. Methods depending on non-linear processing of the gradient along paths in the image fall to smoothly generalize to two dimensions, and cannot predict the appsarance of images in uhich different pethe reeuit In different lightneseee. ON LIGHTNESS 19 2.' Lightness in Two Dimensional [mages. 2.1 Tuo-Oi mansions! flathod - Continous CasSt Ue naed to extend our Ideas to tuo dimensions In order to deal wl th actual images, Thsra are a number of Mays of arriving at tha process to be described here, we shall follow tha simplest (Horn 19681. Urn need to find two-dimensional analogues to differentiation and Integration. The first partial derlvatlvee %r% directional and thus unsuitable alnce they Mill for example completely eliminate evidence of edges running In a direction parallel to their direction of differentiation. Exploring the partial derivatives and their linear combinations one finds that the Laplacian operator la tha louest order combination that la leotroplc, or rotational ly symmetric. Tha Laplacian operator is of course the sum of the second partial derivatives. 2.1.1 Applying the Laplacian to a flondrian. Before investigating the invert ibi I i ty of thie operator* let ue eee yhat happene when one appliee It to the image of a Mondrlan. Inelde any region one Mill obtain a finite value due to the variation In Illumination Intenelty. At each edge one Hi 1 1 get a pulae pair, one positive and one negative. The area of each pulse ulll be equal to the Intenelty step. . -* << 1 > 2 j FIGURE 3: Applying the laplacian operator to the image of a Hondrian figure. ON LIGHTNESS 21 Thia can beat be esen by considsring tha first derivative of a step, nanely a single pulee. If thie ia differentiated again one obtains a doubled pulse as described. Since this pulse Mill extend along the edge, one mag think of it as a pulse-wail. So each edge separating regions m1 1 1 produce a doubled pulse nail. It la clear that one can once again separata the coaponant due to reflectance and i I I umi nation simply by discarding all finite parts. 2.1*2 Inverse of the Laplacian Opsrator. To co«p)ete the task at hand one than has to find a process for undoing the effect of applying the Laplacian. Again there are a number of approaches to this problem, we will use the ahortest IHorn 1968). In essence one has to solve for p(x,y) in a partial differential equation of the formi L(p(x,y)) - d(x,y) Thie ia Poiason'e equation and It la usually solved inside a bounded region using Green's function (Garabsdien 13541 1 «K,y) - MC{L^« t g)td^^)d{di, ON LIGHTNESS 22 The form of Gnm'i function G» dependa on the ehape of the region boundary Nom If the retina ie infinite all polnta are treated similarly and Green* a function dapenda only on tuo para we tare, (J; - x) and (*) - y) . Thla poeit tonal indepandanca lapllee that the above integral almply beconea a convolution. It can ba ehown that Green'e function for thi a caaa lai G(£ t *J * *,y) - Q/2 * ) log f (l/r) Uhere r x - <£,- x) + (*|- y) So p(x ,y) - / / (1/2") log ( (l/r) * dt^,*| ) d^ d^ Thua the Invaraa of the Laplacian operatora la almply convolution ulth (1/2*) log <l/r). To ba preciee one haet (|-=+ i*£ // d/2"> log c (l/r) * d(^,T ) d^dij - d(x,y) Thie ie the two-dimeneional analogue oft ON LIGHTNESS 23 *t fit) dt - f(x) 2*1.3 Why one can use the Convolutional Inverse. If the retina le coneidsred infinite one can express the Inverse ae a elmple convolution. If the retina le finite on the ether hand one hae to uee the acre complicated Creen*8 function foraulation. Now conelder e scene on a uniform background whose Image le totally contained on the retina. The result of applying the forward transform and thresholding will be zero In the area of the uniform background. The convolutional Inverse will therefore receive no contribution from outside the retina. As a result one can use the convolutional form of the Inverse provided the Image of the ecene le totally contained within the ratine. 2.1.4 Normal Ization. Once again one flnde that the reconstructed reflectance le not unlqua. That le. any non-singular solution of Llpfa.y)) - 8 can be added to tha Input without affecting the result. On the Infinite plane auch eolutions have the form p(x t y) » (a*x + b*y + c). If the scene only occupies a finite region of space it can be further shown that the eolution will be unique up to a constant and that one doee not have to ON LIGHTNESS 24 worry about possible slopes. To be specifics the background around the ecene Ml I I be constant in the reconstruction. So one has here exactly the same normalization problem as In the ons-dleenslonal case. Assigning uhlts to the region with highest numerical value In the reconstructed output appears to be e reasonable method. 2*1.6 Tuo-Dimeneional Continous (lethodi Details* Let r'(x.y) be the reflectance of the object at the point corresponding to the iaags point <x t y). Let e* (x,y) be the source intensity at that object point. Let pMx.y) be their product, that ie the intensity at tht iaage point U.y). Note that rMx.y) and a*(x,y) are poei tlve. p* (x,y) - a* (x t y) ft p 1 (x,y) Let p(x a y) be the logarithm of p* (x.y) and so on. p(*.y> - a(x,y) t r(x«y) Now assume that s(x,y) and Its first partial derivatives are continoue - a reasonable assumption to make for the distribution of illumination on the object. Let L be the Laplacian opsrator* ON LIGHTNESS 25 d(x,y) - L(p(x,y)) - Llsfx.y)) + L(r(x,y)) Nom L(i(k»u)I Mill be flnl te everywhere, whl le L(r (x,y) ) mIII bo zero except at each edge eeparating regions, where one nil I flr>d a double pulse wall ae deecribed. Now discard all finite partet t<x,y) - T(L(p<x.y))> -Ltrlx.yM Lot G bo the operator corresponding to convolution by (1/2* ) log,(l/r] Note that (G)" 1 - L. ■U,y) - G(T(L(p(* t y))) - r(x t y) +C 2.2 Tuo-Dlmenslonel flethod - Oiecrete Case. Once again we turn fron a contlnous image to one sampled at dlecrata points. First Me wit l have to decide on a teaselatlon of the Image plane* 2.2.1 Teaselatlon of the .nage Plans. For regular teseelations the choice la between triangular, aquara and hexagonal unit eel la* In much past work on image proceeding, aquar teaeelations have been used for the obvious reaeone. Thle particular taeaalation of the Image has a number of disadvantages. Each cell hae ON LIGHTNESS 26 two kindt* of neighbors, four adjoining the sides, four on the corners. This results In a number of aaymmstr lee. It makes It difficult to find convenient difference schemes approximating the Laplacian operator with low error term. Triangular unit cells are evsn worse In that they have three kinde of neighbors* compounded these drawbacks. Note also that near-circular objecte pack tightsst in a pattern with hexagonal cella. For these reasons hs mIM use a hexagonal unit cell. It should be kept In Bind however that It is easy to develop equivalent results using different tesselat Ions. 2.2.2 Discrats Analogue of tha Laplacian. Having decidsd on ths teaselation we nssd now to find a discrete analogus of the Laplacian operator. Convolution Mlth a central positive value and a rotational ly syuetric negative surround of equal weight le one possibility. Aside from a negative scale factor, this Mill approach application of the Laplacian in the Hull as ths cell size tends to zero* If one were to use complicated surrounds, the trade-offs between accuracy and resolution would suggest using a negative surround that decreasss rapidly outward. For the sake of simplicity we will choose convolution with a central cell of weight 1, surrounded by elx eel le of weight -1/6. Thle function is convenient, symmstrlc end has a small error ten* It Is equsl to - ( h/4 L + h/64 L 1 ) plus sixth and higher FIGURE 4a: A discrete analogue of the laplacian operator. . -*■ Nflwvtti Delta function minus this discrete analogue. ON LIGHTNESS 28 order derivative* IRlchtmeyer 19571 . It should again be pointed out that similar results can be developed for different functions. 2-2.3 Inverse of the Oiscrete Opsrator, The forward differencing operator has the form «»i ■ p., -J}^h * **\ Uhere p,« 1* the logarithm of linage intensity, m ( ; fti'e Might*, uhleh In our case are 1/G, and the sun is taken over the six immediate neighbors* Ue now have to determine the inverse operation that recovers d-- from d|i* One approach ie to try and solve the difference equation of the form v5j Vi ^ p w " -| i Or In matrix form U p • d. Note that U la tparse, hevlng l*e on the diagonal and -1/G'e scattered around. For a finite retina with n eeneor cells one has to introduce boundary conditions to ensure that one has as many equations ae there are unknowns. One then eimply inverte the ON LIGHTNESS 29 matrix U and getst p - U~ d. This is entirely analogous to the solution in the continous case for a finite retina. U~ corresponds to the Green's function. Much as Green's function has a large "support", that is, is non-zero over a large area, so 1*1 is not sparse. This implies that a Jot of computation is needed to perform the inverse operation. 2.2.4 Computational Effort and Simplification. Solving the difference equations for a given image by simple Gauss- Jordan elimination requires of the order of n /2 arithmetic operations. Another approach is to invert U once and for all for a given retina. For each image then one needs only about n x arithmetic operations. Note that the other operations, such as forward differencing, require only about Gftn ar i thmet ic operations. Uhat in effect is happening is that each point in the output depends on each point in the differenced image. Both have n points, so n 1 operations are involved. Not only does one have to do a lot of computation, but must also store up the matrix U of size n . This le quite prohibitive for even a small retina. This latter problem can be avoided if one remembers the simplification attendant to the use of an infinite retina in the continous case. There ne found that the integral ulth Green's function simplified into a convolution. Similarly, if one assumes an infinite retina here, one finds that W and its inverse become very regular, Tho ON LIGHTNESS 38 rows in U are then all the sane and the same is true of U . Each value in the output then depends in the same way on the neighboring points In the differenced image. One need only store up the dependence of one point on its neighbors for this simple convolutional operation. The only remaining difficulty is that U 19 now Infinite and one can no longer invert it numerically - one has to find an analytical expression for the inverse. I have not been able to find this inverse exactly. A good first approximation is log,lr D /rl - except for r » 0, when one uses 1 + log,(r c K Here r is the distance from the origin and p^ is arbitrary. The reaainder left over i-ihen one applies the forward difference scheme to this approximation lies between log.ll + r ) and log 11 - p" ) t This error term is of the order of r • In practice one does not have an infinite retina, but as has been explained for the continous case one can use the convolutional method described above for a finite retina, provided that the image of the scene is wholly contained inside the boundaries of the retina* It is possible to find an accurate inverse of this kind valid for a limited retinal size by numerical means. 2.2.5 Two-Dimensional Discrete Method: Details. Let r*. be the reflectance at the object point corresponding to the image point (l.j). Let »y be the intensity of the incident light at this object point. Let &l be the intensity in the image at point ON LIGHTNESS 31 »?j ■•b-4 Let p;j be the logarithm of pV and so on. Ltt L* be the operator that corresponds to convolution with tha analogua of tha Laplaclan* Lat G* be i ta Inverse. df i * p 'i "2-?^* p w ,d - u( p» J Tha weigthe m- are 1/6 in thte casa t and tha eua la taken over the ah Immediate nelghbora. else Here the sum extends over the whole retina and v-» is the convolut ional Inveree found numerically as explalnsd above. ON LIGHTNESS 32 2.2.G Simplicity of the Invert*. The forward transform. Involving only a slmpls subtraction of Immediate neighbors, Is clearly a rapid, local operation. The inverse on the other hand ie global, since each point in the output depends on each point in the differenced image. Computationally this makes the inverse slow. The inveree is simple in one sense houeveri The difference equations bslng solvsd by ths Inverse have the seme form as the equations used for the forward transform end are thus local. The problem Ie that the output here feeds back into the system and effecte can propagate across ths retina. The apparent global nature of the inverse is thus of d rather special kind find, as we will see later, gives rise to very simple implementations involving only local connect tone. 2.2.7 Iterative Methods of Solution. There are of course other methods for solving large sete of equations. The fact that U is sparse and has large diagonal elements* euggeets trying something I ike Gauss-Ssldsl iteretion. Each Iteration takes about S*n arithmetic operations. For effects to propagate acroee the retina one requires at least ^(4ftn* - l)/3 iteratlone. Thle le becauee a hexagonal retina of width m has (3*m x + l)/4 cells. The above euggeets that one might be able to get away with less than n l arithmetic operatione. In practice It Is found that effecte propagate very elouly ON LIGHTNESS 33 and many more Iterations art needed to stabilize the solution, Ona does not hava to atora U t alnca it la aaailg generated aa one goea along. Iterative echemee correspond to adding a tims-der ivative to tha Poiaaon equation and so turning It into a heat-equation, Aa ona continues to iterate the steady-state solution la approached. This intuitive modal givea some insight into how the process will converge* 2.2*8 Convergence of Iterative and Feed-back Schemes. It is not immediately clear that iterative achamea of solving the difference equations hIJI converge. If they do, they will converge to the correct solution. Let 5 be the delta function, that is. one at tha origin, zero elsewhere. It can be shown that if tha forward convolutional operator la u* the convergence of Iterative echemee depends on the behaviour of the error term, ( 6 * w) , aa n bacomaa large. Raiaing a convolutional operator to an integer power is intended to elgnlfy convolution with itself. In our case, m is ons at the origin, ulth elx valuea of -1/6 around it. So ( 6 - y) ill 1 1 b« zero at the origin with elx valuaa of 1/6 around it. New while ( i- w) Ml 1 1 aluays have a total area of one. It doea epread out and I te value tends to zero at every point aa n tende to Infinity. So thie iterative scheme convsrgeet similar results could bo derived for other negative surrounds. ON LIGHTNESS 34 2.2*9 Setting the Threshold. In the dlecrete case a finite threshold must be selected. As before, let g' be the seal lest stsp In the logarithm of reflectance In the ecene, h the sensor spacing and u the radius of the point-spread function of the optical system* Than us define the effective minimum step asi g • g' * mind, h/2u) There are some minor differences In what follows depending on whether one considers the sensor outputs to bo intensity aanplee at cell-centore or averages over the cell area. The smallest output dus to Mn edge ui II be about g/6. This Is producsd uhsn ths edge le oriented to cover Just one cell of the neighborhood of six. Let Q be the maximum of the intensity gradient - that is the Laplacian of intensity In this case. Choose the threshold s as follows: e < g/G e > * h* e > 3^7/G" * ON UGHTNESS 35 2*2*10 Sons Notes on This tlethod* Notics that an illumination gradient that varies aa soms power of distance across ths image becomes a linear slops after taking logarithms and thus produces no component after the differencing operations* Such simple gradients are eupprseeed even without the thresholding operation. In practice ths parameters used in choosing the threshold may not be known or may be variable. In this cass ons can look at a histogram of the differenced image. It will contain values both positive and negative corresponding to edges and also a large number of values clustered around zsro due to i I lunination gradients, noise and so on* The threshold can bs convsnisntlu choosen to contain this central blob* Noise and illumination gradients havs an effect similar to that In the one-dimensional cass. Ulth finite cell spacing ons cannot precisely eeparate ths two components of ths image intensity and at each edge the information Mill be corrupted slightly by nolss and Illumination gradient. As the density of sdgas per csll arsa goss up the effect of this becomes mors apparent. In highly textured scenes the illumination gradient is hard to eliminate* Once again one hee to dsclde on a normalization echeme. The best method probably Is to 1st the highest numerical value In the ■ reconstructed output correspond to unite* ON LIGHTNESS 36 2.2.11 Dynamic Range Reduction. Applying the retinex operation to an image coneiderable reducee the range of vatuee. Th'io ie because the output, being related to reflectance, will only have a range of one to two ordere of magnitude, while the Input will a I eo have i lluminat ion gradients. This will make euch processing useful for picture recording and transmission (Horn 19681* 2.2.12 A Frequency Domain Interpretation* t It may be of Intoreftt to look at thit as t hod fron got another point of view. Uhat one does Is to exentuate the high-frequency components, threshold and then attenuate the high-frequency coeponente. To eee this, coneider first the forward operation. The fourier traneform of the convolut tonal operator corresponding to differentiation le 1«. Similarly the two-dimensional fourier traneform of the convolutional operator corresponding to the LapJeclan le - p*. Here p ie the radiue in a polar coordinate system of the two-dimensional frequency space. In either case one is multiplying the fourier transform by sone function th^t increaeee with frequency. Now consider the revsree operation. The fourier traneform of the convolutional operation corresponding to Integration ie l/l«. Similarly the fourier transform of (1/2*) log r (1/r) ie -l/p x . So In the Inverse step one undoes exactly the emphaele given to high frequency components in the forward operation. ON LIGHTNESS 37 In both thn ons-dlmenslonal and the tuo-dinensional case one loeee the zero frequency component. Thle It why the reeult hae to be normal ized. 2.3 Phgelcal llodele. There are nueerous continous physical nodele to illustrate the inverse traneformatlon. Anything that satisf lee Polsson'a equation ulll do. Such physical models help one visualize Mhat the Inverse of a given function might be. Examples in two dimensions arei perfect fluld-f Iom» eteady diffusion, steady heat-flou, deformation of an elastic membrane, electro-statics and current flow in a resistive sheet. In the last model for example, the input is the distribution of current flowing Into the reeietlve sheet noraal to Its surface, the output ie the dletrlbutlon of electrical potential over the surface. In addition to helping one visualize solutions* these contlnoue eodels aleo suggest discrete models. These can be arrived at elaply by cutting up the two-dimensional space in a pattern corresponding to the interconnection of neighboring cells. That Is, the remaining parte form a pattern dual to that of the sensor cell pettern. Us ulll discuss only one such discrete mods It ON LIGHTNESS 38 2.3.1 A Discrete Physical node I. Consider the reeletive eheet described* cut up in the dual pattern of the hexagonal unit cell pattern. What will be left ie an Interconnection of resistors In a triangular pattern. The inpute to thie system will be currents injected at the nodes, the potential at the nodes being the output. Thie then provides a very simple analog implementation of the tedious inverse computation. It ie perhaps at first surprising to sss that each cell ie not connected to every other In a direct fashion. One would expect thie from the form of the computational Inverse. Each cell In the output does Of COUrft* h«V* ft connection Via the ftthtr ceils to oich of thi Inputs. Paths are shared however in a way that makes the result both simple and planar. Conelder for the moment just one node. The potential at the node ie the average of the potential of the six nodes connected to it plue the current injected times H/B, where R Is the resistance of each resistor. The economy of connection is due to the fact that the outpute of thie system are fed back into it. It also illuetrates that thie model locally solves exactly the same difference equation ae that used In the forward transform, only now in reverse. This immediately suggests an important property of thie model! By eimply changing the Interconnections one can make an inverse for other forward transforms. Slmpleet of all are other image plane tesse lat ions, both regular and irregular. One simply connecte the resistors in the ■ '* FIGURE 5: Resistive model of the inverse computation. The inputs are the currents injected at the nodes. The outputs are the potentials at the nodes* ON LIGHTNESS 40 i sane pattern as are tha calls in tha input* Mora complicated weighted surrounds can ba handled bg using resistors with resistances inversely proportional to tha weights* Tha network of resistors will then no longer be planar. 2*3.2 A Feed-back Scheme for the Inveree. Both the comment about outputs feeding back into the resistive Model and the earlier notes about iterative schemes suggest get another interesting model for the inverss using linear cunning devices. Operational amplifiers can serve this purpose. One simply connects the summing element so that they solve the difference equation implied by the forward transform. Once again It la clear that such a scheme can be generalized to arbitrary tesselations and wsighted negative eurrounde eimply by changing the Interconnections and attenuations on each input. Some questions of stability arise with esotsrlc Interconnections. For the simple ones stability Is assured. A little thought will show that ths resistive model described earlier ie in fact a more economical implementation of just this scheme with the difference that there the inputs are currents, while here they are potentials. £ u *j id .a TJ c i m "E ■u 2 c t. <0 U1 4-t e 2 i +J a 3 c^ ■ c: M- *C o c £ ^ o VI 40 •*— c *> <v "- ^ *- o *-* ft) 8 V ^ ^ ^ -j 1- CJ I > *^ ON LIGHTNESS 42 2.4 Limitations of the Simple Scheme Presented. Tha Method presented hara ul II not correctly calculate reflectance if uaed unmodified on general scenes. It may however calculate lightness fairly well. Ae the msthod stands now for example, a sharp ehadow edge ulll not be distinguished from a real edge In the scene and the two regions so formed nil 1 produce different outputs, whi le their reflectances are the same. It nay be that thie le reasonable nevertheleee, elnce ue perceive a difference in apparent lightness. Smooth gradations of reflectance on a surface due either to shading or variations In surface reflectance ulll be eliminated by the thresholding operations except as far as they affect the intensity at tha borders of the region* Thie may imply that ue need additional channele In our visual eystem to complemsnt the onee cerrylng the retlnened Information elnce ue do utilize shading as a depth-cue* The simple normalization scheae described Mill also be sensitive to specular reflections, flourescent paints and light-sources in the field of view* Large depth-discontinul ties present another problem. One cannot assume that the illumination is equel on both sides of the obscuring edge. In this case tha illuminating component does not vary smoothly over the retina, having inetead some sharp edges. S -2 C ** *-> 3 £ E as .woe E ^ 0> r- 5 2 '2 i E 5 4 5 V a > in o Id CI a. o L/l 4B 3 It 1% o L3 ON LIGHTNESS 44 2.S Computer Simulation of the Discrete Method* A computer program was used to simulate the retinex process described on a small retina with both artificial and real images eeen through an image dissector camera. The hexagonal unit cell le used in this program and the retina itself Is also hexagonal. The retina containe 1027 cells in a pattern 36 celle across. Thie ie a compromise dictated by the need to limit the number of arithmetic operatlone In tha inveree traneform. In this case one needs about a million and thie takee about a minute of central processor time on our PDP-18. Both the art 1 f Iclal and the real Hondrians consist of regions bounded by curved outline* to smphaaiza that this ntthod does not require straight-line edges or boundsry extraction and description. Various distributions of Incident illumination can be selected for the artificial ecenes. In each case the processing satisfactorily removes the gradient. For the real scenes It is hard to produce really large Illumination gradiente by poeltloning the light-sources. The reconstruction does eliminate the gradient well, but often minor flaws m1 1 1 appear In the output due to noiee in the input and a number of probleme ulth this kind of Input device such as a very considerable scatter. It ie not easy to predict uhat effects euch imaging device defects will have. The output is displayed on a OEC 349 display which has a mere eight grey-levele. It uould be interesting to experiment ulth larger retinae and better Image Input- and output-devices. ON LIGHTNESS AS 2.5.1 For« of Inverse used In the Computer Simulation* The convo I u t i one I form of the Inverse was used for speed and lou storage requirement. This neccesitated solving the difference equation! once, given a pulse as Input. The symmetry of the hexagonal pattern allous one to identify symmetrically placed cells and only 324 unknowns needed to be found for a convolutional inverse sufficient for the eize of rotina described. As Mentioned before, this function is closely approximated by log, (r # /rl for large r. This can be used to establish boundary condi t ions. FIGURE 8: The method applied to an artificial image. FIGURE 9: The method applied to a real image FIGURE 10: The method applied to Craik's figure. FIGURE 11: Apparent lightness predicted for incomplete figure- The subfigures in the above have the following interpretation: A Input - logarithm of image intensity B Differenced image C Thresholded difference D Output - computed lightness E Illumination distribution (P^*!**) 1 fj ■i.i ON LIGHTNESS 47 3* Implications and Conclusions. 3*1 Parallel Imaga Procsaaing Hardware. Tha methods described hare for forward transforming, thresholding and Inverse transforming inaediately tempt one to think In terms of electronic componenta arranged in parallel layere. Enough has been aald about different models to make it clear how one Bight connect euch components. Large scale integrated circuit technology may be useful* provided the signals are either convsrtsd from analog to digital for* or batter still, good linear circuits are available In thia form. Construction of such devlcea would be premature until further experimentation le performed to decide on optimal tesselatlone, optimal negative surrounds, thresholding operations and normalization schemes. Theee decleions are best guided bg computer simulation. 3.2 Cognitive Peychology. One of the artificial ecanes was created to Illustrate Cralk's lllueion (Brindley I960, Cornsweet 1970). Here a sharp edge le bordered by eecond-order gradients. As one might expect the smooth gradlente are loet in the thresholding and reconstruction produces two regions each of uniform brightness* Ths diffsrencs in brightnsss between the reglone le equal to the original Intensity stsp at the edge. The fact that the process presented hsre felle prey to this ON LIGHTNESS 48 i I lueion is of course no proof that humans use the «a»* mechanism. It le interesting that this technique allows one to predict for example the appearance of pictures containing incompletely closed curves with second-order gradients on either side. 3*3 Neuro-physr ology. The method described here for obtaining lightness from Image Intmnelty euggeete functions for a number of etructuree in the primate retina. The horizontal ceils appear to be involved In the forward transformation, while some of the amacrlnes may be Involved in the inverse transformation* For details see tha paper by David Narr IDarr 19741 , In which he uses this hypothssis to explain an astonishing number of facte about the retina. 3.4 Concluelon. A elmple layered, parallel technique for computing lightness from image Inteneity has been preesntsd. Ths method does not involve an ability to describe or understand the ecene, relying instead on the spatial differences In the distribution of reflectance and illumination. The forward step involves accentuating ths edges between regions. The output of thle step is then thresholded to remove Illumination gradiente and noise. The Inverse step merely undoes the accentuation of the edges* ON LIGHTNESS 49 Physical models nave been glvan uhich can perform this commutation efficiently in parallel layers of simple networks. The method haa been simulated and applied to a number of images* The method grew out of an attempt to extend Land'e method to two dimenslone and fills the need for a llghtneee Judging process in his retinex theory of color perception. The possibility of processing an imags In such a parallel* elmpla fashion uithout higher-level understanding of the scene reinforces my belief that such lou-level processing is of Importance In dealing ulth a number of features of Images. Amongst thsee are shading, stereo disparity, focue, edge detection, scene segmentation and motion parallax. Some of thia kind of processing mag actually happen In the primate retina and viaual cortex. The Implications for Image analysis are that it may well be that a number of such pre-proceeeing operatlona should be performed automatically for the whols Image to accentuate or extract certain attributes before one brings to bear the more powerful* but tedious and slou sequential goal-directed methods. CN LIGHTNESS SB B lb I iographg. Brindleu. G.S, (1368) "Physiology of the Retina and Visual Pathway." (Honograph No. 6 of the Physiological Society). Londoni Edward Arnold Ltd. Cornsueet, T. (1978) "Visual Perception." Nau Vorki Academic Preae. Garabedian, O.R. (19G4) 'Partial Olffarantlal Equations." Nau Yorki John Ui lay. Goetha, J. 14. von. (1818) "Zur Farbenlehre. " Tuebingen. Hardy, A.C. (Ed,) (193S) "Tha Handbook of Coioriaetry. " Cambridge, Haest n. I ,T» Praaa* Hefmhoitz, H.L.F. (1867) "Handbuch dar Physlologlachan Optlk." Lalpzlgi Voaa. Alao tranalatedi Southall, J.P.C.S. "Handbook of Physiological Optics." Nau Vorki Oovar Publications. He I son, H. (1938) "Fundamental Problems In Color Vision I." Journal of Experimental Psychology, 23. Helaon, H. (1948) "Fundamental Problems in Color Vision II." Journal of Experimental Psychology, 26. Herlng, E. (1875) "Zur Lahrs von Lichtsinne - Grunzuege elner Theorie des Farbsinnes." SBK. Akad. Uiss. Uien ftath. Naturuies. K.7B. Horn, B.K.P, (1368) "Tha Application of Fourier Transform Methods to Image Procasalng." S.M. Thesis. E.E. Department, fl.I.T. pp. 81-87 A pp. 93-94. Judd, D.B* (1948) "Hue, Saturation and Lightness of Surface Colore with Chromatic Illumination ." Journal of the Optical Society of America, 38. Judd, D.B. (1962) "Color In Business, Science and Industry." New Vorki John Ulley. Land, E.H. (1959) "Experlmsnte in Color Vieion," Scientific American. Land, E.H, (1964) "The Retinex," American Scientist. S2. Land, E.H. 4 flcCann, J.J. (1971) "Lightness Theory." Journal of the Optical Society, 61. Lettvin, J.Y. (1967) 'The Color of Colored Things." Quarterly Progress Report, 87, Rssearch Laboratory for Electronics, M.I.T. ON LIGHTNESS 51 Marr f D- (1974) "An Analysis of the Primate Retina." A.I. Clemo, 296, Artificial Intel I Igence Laboratory, M.I.T. tlaxueH, J.C. (185G) "On the Unequal Sensitivity of the Foramen Centrale to Light of Different Colours. " Rep. Brit. Assoc. Newton, Sir Isaac. (1784) "OpticKs." Londoni Sanuel Smith 6 Benjamin Uatford. Also Neu Yorkt Dover Publ icatione. Richards, U. (1971) "One-stage flodsl for Color Convereion. ■ Journal of the Optical Society, 62. Rlchtmeyer, R.O. & dor ton, K.U. (1957) "Difference Hethode for Initial Value Problems.- Neu YorKi John U isy. Young, T. (1828) "On the Theory of Light and Colour." Philosophical Transact lone* Aleo im Teevan, R.C. & Birney, R.C. (Ed.) (1361) "Color Vision," Van Noetrand.